import numpy as np
from typing import Dict, List, Tuple
from physics import PhysicsModel, PhaseType
from parallel_compute import (
    calculate_mass_flux_array,
    calculate_momentum_flux_array,
    calculate_wall_friction_array,
    calculate_interfacial_friction_array,
    calculate_lateral_flow_resistance_array,
    calculate_turbulent_mixing_array,
    solve_momentum_equation_parallel,
    solve_continuity_equation_parallel
)

class Solver:
    """两相流求解器类（使用并行计算）"""
    def __init__(self, grid, physics_model):
        self.grid = grid
        self.physics = physics_model
        
        # 求解器参数
        self.max_iterations = 1000
        self.convergence_tolerance = 1e-6
        self.relaxation_factors = {
            'pressure': 0.3,
            'velocity': 0.7,
            'void_fraction': 0.3
        }
        
        # CFL条件相关
        self.max_cfl = 0.8
        self.min_dt = 1e-6
        self.max_dt = 0.01
        
        # 时间步长
        self.dt = self.min_dt
        
        # 设置初始条件
        self.initialize_fields()
        
    def initialize_fields(self):
        """初始化所有场量"""
        print("\n初始化计算场...")
        for subchannel in self.grid.subchannels:
            n_axial = len(subchannel.alpha)
            
            # 1. 压力场初始化：线性分布（常压，略高于大气压）
            P_in = 1.2e5  # 入口压力 (1.2 bar)
            P_out = 1.0e5  # 出口压力 (1.0 bar)
            subchannel.pressure = np.linspace(P_in, P_out, n_axial)
            
            # 2. 空泡份额初始化：入口值延拓
            alpha_in = 0.3  # 入口空泡份额（减小初始空泡份额以提高稳定性）
            subchannel.alpha = np.full(n_axial, alpha_in)
            
            # 3. 速度场初始化（降低初始速度以提高稳定性）
            v_l_in = 0.5  # 入口液相速度
            v_g_in = 0.8  # 入口气相速度（略快于液相）
            subchannel.velocity = {
                'liquid': {
                    'axial': np.full(n_axial, v_l_in),
                    'lateral': np.zeros(n_axial)
                },
                'vapor': {
                    'axial': np.full(n_axial, v_g_in),
                    'lateral': np.zeros(n_axial)
                }
            }
            
            # 4. 保存初始状态
            subchannel.save_current_state()
            
        print("场初始化完成：")
        print(f"- 压力范围：{P_in/1e5:.2f} bar - {P_out/1e5:.2f} bar")
        print(f"- 入口空泡份额：{alpha_in:.3f}")
        print(f"- 入口速度：液相 {v_l_in:.2f} m/s, 气相 {v_g_in:.2f} m/s")
        
    def calculate_time_step(self) -> float:
        """计算自适应时间步长"""
        # 计算最大速度
        max_velocity = 1e-6  # 设置一个最小速度以避免除零
        for subchannel in self.grid.subchannels:
            for phase in ['liquid', 'vapor']:
                max_velocity = max(max_velocity,
                                 np.max(np.abs(subchannel.velocity[phase]['axial'])),
                                 np.max(np.abs(subchannel.velocity[phase]['lateral'])))
        
        # 基于CFL条件计算时间步长
        dt_cfl = self.max_cfl * min(self.grid.dx, self.grid.dz) / max_velocity
        
        # 限制时间步长在允许范围内
        self.dt = np.clip(dt_cfl, self.min_dt, self.max_dt)
        
        return self.dt
    
    def predict_velocity(self, subchannel, dt: float):
        """速度预测步骤"""
        n_axial = len(subchannel.alpha)
        
        for i in range(1, n_axial-1):
            # 获取物性
            props_g = self.physics.calculate_properties(subchannel.pressure[i], PhaseType.VAPOR)
            props_l = self.physics.calculate_properties(subchannel.pressure[i], PhaseType.LIQUID)
            
            # 使用动量方程预测速度
            for phase, props in [('vapor', props_g), ('liquid', props_l)]:
                # 计算压力梯度
                dp_dz = (subchannel.pressure[i+1] - subchannel.pressure[i-1]) / (2 * self.grid.dz)
                
                # 预测速度
                subchannel.velocity[phase]['axial'][i] = (
                    subchannel.previous_velocity[phase]['axial'][i] -
                    dt * dp_dz / props['density']
                )
    
    def solve_pressure_correction(self, subchannel):
        """求解压力修正方程"""
        n_axial = len(subchannel.alpha)
        
        # 构建压力修正方程系数矩阵
        A = np.zeros((n_axial, n_axial))
        b = np.zeros(n_axial)
        
        # 设置入口边界条件（固定压力）
        A[0, 0] = 1.0
        b[0] = 0.0  # 入口压力不需要修正
        
        # 设置出口边界条件（固定压力）
        A[-1, -1] = 1.0
        b[-1] = 0.0  # 出口压力不需要修正
        
        # 构建内部节点的方程
        for i in range(1, n_axial-1):
            # 获取物性
            props_g = self.physics.calculate_properties(subchannel.pressure[i], PhaseType.VAPOR)
            props_l = self.physics.calculate_properties(subchannel.pressure[i], PhaseType.LIQUID)
            
            # 计算密度项（添加保护）
            rho_g = max(props_g['density'], 1e-6)
            rho_l = max(props_l['density'], 1e-6)
            
            # 计算速度项
            v_g = subchannel.velocity['vapor']['axial'][i]
            v_l = subchannel.velocity['liquid']['axial'][i]
            
            # 计算系数（考虑密度和空泡份额的影响，添加保护）
            alpha = np.clip(subchannel.alpha[i], 1e-6, 1-1e-6)
            coeff = (alpha / rho_g + (1 - alpha) / rho_l)
            
            # 添加数值稳定项
            stability_factor = 1e-8
            
            # 设置系数矩阵
            A[i, i-1] = -coeff / (self.grid.dz ** 2)
            A[i, i] = 2 * coeff / (self.grid.dz ** 2) + stability_factor
            A[i, i+1] = -coeff / (self.grid.dz ** 2)
            
            # 计算质量不平衡作为源项
            G_g = alpha * rho_g * v_g
            G_l = (1 - alpha) * rho_l * v_l
            
            # 计算质量流量的空间导数（使用中心差分）
            alpha_prev = np.clip(subchannel.alpha[i-1], 1e-6, 1-1e-6)
            v_g_prev = subchannel.velocity['vapor']['axial'][i-1]
            v_l_prev = subchannel.velocity['liquid']['axial'][i-1]
            
            G_g_prev = alpha_prev * rho_g * v_g_prev
            G_l_prev = (1 - alpha_prev) * rho_l * v_l_prev
            
            dG_g_dz = (G_g - G_g_prev) / self.grid.dz
            dG_l_dz = (G_l - G_l_prev) / self.grid.dz
            
            b[i] = -(dG_g_dz + dG_l_dz)
        
        try:
            # 添加对角优势以提高稳定性
            for i in range(n_axial):
                A[i,i] = max(A[i,i], 1.1 * np.sum(np.abs(A[i,:])))
            
            # 使用直接求解器求解压力修正方程
            pressure_correction = np.linalg.solve(A, b)
            
            # 检查解的有效性
            if not np.all(np.isfinite(pressure_correction)):
                raise ValueError("压力修正解包含无效值")
            
            # 限制压力修正的幅度
            max_correction = 0.1 * np.mean(subchannel.pressure)
            pressure_correction = np.clip(pressure_correction, -max_correction, max_correction)
            
            # 应用压力修正（带松弛因子）
            subchannel.pressure += self.relaxation_factors['pressure'] * pressure_correction
            
            # 确保压力始终在物理范围内
            subchannel.pressure = np.maximum(subchannel.pressure, 1e5)  # 最小压力为1个大气压
            
        except (np.linalg.LinAlgError, ValueError) as e:
            print(f"警告：压力修正方程求解失败: {str(e)}")
            # 使用简单的压力梯度估计
            dp = (subchannel.pressure[-1] - subchannel.pressure[0]) / (n_axial - 1)
            subchannel.pressure = np.linspace(subchannel.pressure[0], subchannel.pressure[-1], n_axial)
    
    def correct_velocity(self, subchannel):
        """基于压力修正更新速度场"""
        n_axial = len(subchannel.alpha)
        
        for i in range(1, n_axial-1):
            for phase in ['vapor', 'liquid']:
                # 计算速度修正
                dp_correction = (subchannel.pressure[i+1] - subchannel.pressure[i-1]) / (2 * self.grid.dz)
                velocity_correction = -dp_correction * self.relaxation_factors['velocity']
                
                # 更新速度
                subchannel.velocity[phase]['axial'][i] += velocity_correction
    
    def solve_subchannel_parallel(self, subchannel, dt: float):
        """并行求解单个子通道，包含横向流动和空隙漂移"""
        n_axial = len(subchannel.alpha)
        
        # 获取相邻子通道信息
        adjacent_channels = self.grid.get_adjacent_channels(subchannel.id)
        
        # 获取几何参数
        geom = subchannel.get_geometry_parameters()
        
        # 准备数组数据
        alpha = subchannel.alpha
        pressure = subchannel.pressure
        v_g_axial = subchannel.velocity['vapor']['axial']
        v_l_axial = subchannel.velocity['liquid']['axial']
        v_g_lateral = subchannel.velocity['vapor']['lateral']
        v_l_lateral = subchannel.velocity['liquid']['lateral']
        
        # 获取物性
        props_g = self.physics.calculate_properties(pressure[0], PhaseType.VAPOR)
        props_l = self.physics.calculate_properties(pressure[0], PhaseType.LIQUID)
        rho_g = props_g['density']
        rho_l = props_l['density']
        
        # 计算空隙漂移速度
        v_gj, C0 = self.physics.calculate_drift_velocity(rho_g, rho_l)
        
        # 修正气相速度（考虑漂移速度）
        j = alpha * v_g_axial + (1 - alpha) * v_l_axial  # 混合体积通量
        v_g_axial = C0 * j + v_gj
        
        # 轴向动量方程
        for i in range(1, n_axial-1):
            # 计算轴向压力梯度
            dp_dz = (pressure[i+1] - pressure[i-1]) / (2 * self.grid.dz)
            
            # 计算轴向壁面摩擦力
            F_w_g = self.physics.calculate_wall_friction(
                PhaseType.VAPOR, rho_g, v_g_axial[i],
                geom['hydraulic_diameter'], props_g['viscosity']
            )
            F_w_l = self.physics.calculate_wall_friction(
                PhaseType.LIQUID, rho_l, v_l_axial[i],
                geom['hydraulic_diameter'], props_l['viscosity']
            )
            
            # 计算界面摩擦力
            F_i = self.physics.calculate_interfacial_friction(
                alpha[i], rho_g, rho_l, v_g_axial[i], v_l_axial[i],
                geom['hydraulic_diameter']
            )
            
            # 横向动量交换（遍历相邻通道）
            M_lat_g_total = 0
            M_lat_l_total = 0
            for adj_ch in adjacent_channels:
                # 获取相邻通道数据
                v_g_adj = adj_ch.velocity['vapor']['axial'][i]
                v_l_adj = adj_ch.velocity['liquid']['axial'][i]
                
                # 获取连接参数
                conn_params = subchannel.get_connection_parameters(adj_ch)
                
                # 计算横向动量交换
                M_lat_g, M_lat_l = self.physics.calculate_lateral_momentum_exchange(
                    alpha[i], rho_g, rho_l,
                    v_g_axial[i], v_l_axial[i],
                    v_g_adj, v_l_adj,
                    conn_params['gap_length'],
                    conn_params['gap_width']
                )
                
                # 计算湍流交混
                w_t = self.physics.calculate_turbulent_mixing(
                    self.physics.BETA,
                    conn_params['gap_width'],
                    rho_g, rho_l,
                    v_g_axial[i], v_l_axial[i]
                )
                
                M_t = self.physics.calculate_turbulent_mixing_momentum(
                    w_t, v_g_axial[i], v_g_adj
                )
                
                M_lat_g_total += M_lat_g + M_t
                M_lat_l_total += M_lat_l + M_t
            
            # 更新轴向速度
            v_g_axial[i] += dt * (
                -dp_dz/rho_g + F_w_g + F_i/alpha[i] + 
                M_lat_g_total/(rho_g * geom['area']) - self.physics.g
            )
            
            v_l_axial[i] += dt * (
                -dp_dz/rho_l + F_w_l - F_i/(1-alpha[i]) + 
                M_lat_l_total/(rho_l * geom['area']) - self.physics.g
            )
        
        # 更新横向速度（使用横向压力梯度）
        for adj_ch in adjacent_channels:
            conn_params = subchannel.get_connection_parameters(adj_ch)
            dp_dx = (adj_ch.pressure - pressure) / conn_params['centroid_distance']
            v_g_lateral += dt * (-dp_dx/rho_g)
            v_l_lateral += dt * (-dp_dx/rho_l)
        
        # 更新空泡份额（考虑漂移速度的影响）
        j_drift = alpha * (1 - alpha) * v_gj
        alpha_new = alpha + dt * (
            -np.gradient(alpha * v_g_axial, self.grid.dz) +
            np.gradient(j_drift, self.grid.dz)
        )
        
        # 限制空泡份额在物理范围内
        alpha_new = np.clip(alpha_new, 0.0, 1.0)
        
        # 更新结果
        subchannel.alpha = alpha_new
        subchannel.velocity['vapor']['axial'] = v_g_axial
        subchannel.velocity['liquid']['axial'] = v_l_axial
        subchannel.velocity['vapor']['lateral'] = v_g_lateral
        subchannel.velocity['liquid']['lateral'] = v_l_lateral
        
        # 计算残差
        residuals = {
            'velocity_g': np.max(np.abs(v_g_axial - subchannel.previous_velocity['vapor']['axial'])),
            'velocity_l': np.max(np.abs(v_l_axial - subchannel.previous_velocity['liquid']['axial'])),
            'alpha': np.max(np.abs(alpha_new - subchannel.previous_alpha))
        }
        
        return residuals
    
    def solve_time_step(self):
        """求解一个时间步"""
        # 计算时间步长
        dt = self.calculate_time_step()
        print(f"\n时间步长: {dt:.3e} s")
        
        # 预计算物性
        P_min = min(ch.pressure.min() for ch in self.grid.subchannels)
        P_max = max(ch.pressure.max() for ch in self.grid.subchannels)
        print(f"压力范围: {P_min/1e6:.3f} MPa - {P_max/1e6:.3f} MPa")
        self.physics.precalculate_properties(P_min, P_max, n_points=10)
        
        # 保存当前状态
        for subchannel in self.grid.subchannels:
            subchannel.save_current_state()
        
        # SIMPLE算法主循环
        print("\nSIMPLE算法迭代:")
        for iter in range(self.max_iterations):
            max_residuals = {
                'velocity_g': 0.0,
                'velocity_l': 0.0,
                'alpha': 0.0,
                'pressure': 0.0
            }
            
            # 对每个子通道并行求解
            for subchannel in self.grid.subchannels:
                # 1. 速度预测
                self.predict_velocity(subchannel, dt)
                
                # 2. 求解压力修正方程
                self.solve_pressure_correction(subchannel)
                
                # 3. 修正速度场
                self.correct_velocity(subchannel)
                
                # 4. 并行求解守恒方程
                residuals = self.solve_subchannel_parallel(subchannel, dt)
                
                # 更新最大残差
                for key in residuals:
                    max_residuals[key] = max(max_residuals[key], residuals[key])
            
            # 每10次迭代输出一次残差
            if iter % 10 == 0:
                print(f"迭代 {iter:4d}: " + 
                      f"气相速度残差={max_residuals['velocity_g']:.3e}, " +
                      f"液相速度残差={max_residuals['velocity_l']:.3e}, " +
                      f"空泡份额残差={max_residuals['alpha']:.3e}")
            
            # 检查收敛性
            if all(r < self.convergence_tolerance for r in max_residuals.values()):
                print(f"\n已收敛，总迭代次数: {iter + 1}")
                break
            elif iter == self.max_iterations - 1:
                print(f"\n警告：达到最大迭代次数 {self.max_iterations}，未收敛")
        
        # 输出当前状态
        center_channel = self.grid.subchannels[0]  # 使用第一个通道作为示例
        print(f"\n当前状态:")
        print(f"入口压力: {center_channel.pressure[0]/1e6:.3f} MPa")
        print(f"出口压力: {center_channel.pressure[-1]/1e6:.3f} MPa")
        print(f"入口空泡份额: {center_channel.alpha[0]:.3f}")
        print(f"入口气相速度: {center_channel.velocity['vapor']['axial'][0]:.3f} m/s")
        print(f"入口液相速度: {center_channel.velocity['liquid']['axial'][0]:.3f} m/s")
        
        # 清空物性缓存，准备下一个时间步
        self.physics.clear_cache()
        
        # 更新时间
        self.grid.advance_time(dt) 
    
    def solve_piso(self, subchannel, dt):
        """使用PISO算法求解一个时间步
        Args:
            subchannel: 子通道对象
            dt: 时间步长
        """
        # 第一步：预测器步骤
        self.predict_velocity(subchannel, dt)
        
        # 第二步：第一次压力修正
        self.solve_pressure_correction(subchannel)
        self.correct_velocity(subchannel)
        
        # 第三步：第二次压力修正（考虑显式修正项）
        for i in range(1, len(subchannel.alpha)-1):
            # 获取物性
            props_g = self.physics.calculate_properties(subchannel.pressure[i], PhaseType.VAPOR)
            props_l = self.physics.calculate_properties(subchannel.pressure[i], PhaseType.LIQUID)
            
            # 计算显式修正项
            H_g = self.calculate_explicit_correction(subchannel, i, 'vapor', props_g)
            H_l = self.calculate_explicit_correction(subchannel, i, 'liquid', props_l)
            
            # 更新速度
            subchannel.velocity['vapor']['axial'][i] += H_g * dt
            subchannel.velocity['liquid']['axial'][i] += H_l * dt
        
        # 第四步：最终压力修正
        self.solve_pressure_correction(subchannel)
        self.correct_velocity(subchannel)
        
        # 更新空泡份额
        self.update_void_fraction(subchannel, dt)
    
    def calculate_explicit_correction(self, subchannel, i, phase, props):
        """计算显式修正项
        Args:
            subchannel: 子通道对象
            i: 节点索引
            phase: 相态
            props: 物性值
        Returns:
            H: 显式修正项
        """
        # 获取几何参数
        geom = subchannel.get_geometry_parameters()
        
        # 计算对流项
        v = subchannel.velocity[phase]['axial'][i]
        dv_dz = (subchannel.velocity[phase]['axial'][i+1] - 
                 subchannel.velocity[phase]['axial'][i-1]) / (2 * self.grid.dz)
        convection = -v * dv_dz
        
        # 计算扩散项
        d2v_dz2 = (subchannel.velocity[phase]['axial'][i+1] - 
                   2*subchannel.velocity[phase]['axial'][i] + 
                   subchannel.velocity[phase]['axial'][i-1]) / (self.grid.dz**2)
        diffusion = props['viscosity'] / props['density'] * d2v_dz2
        
        # 计算源项（包括重力）
        source = -self.physics.g
        
        return convection + diffusion + source
    
    def update_void_fraction(self, subchannel, dt):
        """更新空泡份额
        Args:
            subchannel: 子通道对象
            dt: 时间步长
        """
        n = len(subchannel.alpha)
        alpha_new = np.zeros(n)
        
        # 获取物性
        props_g = self.physics.calculate_properties(subchannel.pressure[0], PhaseType.VAPOR)
        props_l = self.physics.calculate_properties(subchannel.pressure[0], PhaseType.LIQUID)
        
        # 计算漂移速度
        v_gj, C0 = self.physics.calculate_drift_velocity(props_g['density'], 
                                                       props_l['density'])
        
        for i in range(1, n-1):
            # 计算混合速度
            j = (subchannel.alpha[i] * subchannel.velocity['vapor']['axial'][i] + 
                 (1-subchannel.alpha[i]) * subchannel.velocity['liquid']['axial'][i])
            
            # 计算漂移通量
            j_drift = subchannel.alpha[i] * (1-subchannel.alpha[i]) * v_gj
            
            # 更新空泡份额
            alpha_new[i] = subchannel.alpha[i] + dt * (
                -np.gradient(subchannel.alpha[i] * subchannel.velocity['vapor']['axial'][i], 
                           self.grid.dz) +
                np.gradient(j_drift, self.grid.dz)
            )
        
        # 边界条件
        alpha_new[0] = subchannel.alpha[0]
        alpha_new[-1] = alpha_new[-2]
        
        # 限制空泡份额在物理范围内
        subchannel.alpha = np.clip(alpha_new, 0.0, 1.0) 
    
    def solve_improved_piso(self, subchannel, dt):
        """改进的PISO算法求解一个时间步
        包含：
        1. 动量预测
        2. 压力修正（多次）
        3. 显式修正
        4. 空泡份额修正
        """
        # 1. 动量预测步骤（使用Adams-Bashforth格式）
        if hasattr(subchannel, 'previous_acceleration'):
            # 二阶精度时间推进
            for phase in ['vapor', 'liquid']:
                v_old = subchannel.previous_velocity[phase]['axial']
                a_old = subchannel.previous_acceleration[phase]
                a_current = self.calculate_acceleration(subchannel, phase)
                
                subchannel.velocity[phase]['axial'] += dt * (
                    1.5 * a_current - 0.5 * a_old
                )
                
                # 保存当前加速度
                subchannel.previous_acceleration[phase] = a_current
        else:
            # 首次计算使用一阶格式
            self.predict_velocity(subchannel, dt)
            
            # 初始化加速度存储
            subchannel.previous_acceleration = {
                'vapor': np.zeros_like(subchannel.alpha),
                'liquid': np.zeros_like(subchannel.alpha)
            }
        
        # 2. 多次压力修正循环
        max_pressure_iterations = 3
        for p_iter in range(max_pressure_iterations):
            # 2.1 构建压力修正方程（使用Rhie-Chow插值）
            pressure_correction = self.solve_pressure_correction_rhie_chow(subchannel)
            
            # 2.2 更新压力场
            subchannel.pressure += self.relaxation_factors['pressure'] * pressure_correction
            
            # 2.3 速度修正
            self.correct_velocity_with_rhie_chow(subchannel, pressure_correction)
            
            # 2.4 检查收敛性
            if np.max(np.abs(pressure_correction)) < self.convergence_tolerance:
                break
        
        # 3. 显式修正步骤
        self.apply_explicit_corrections(subchannel, dt)
        
        # 4. 空泡份额修正（使用TVD格式）
        self.update_void_fraction_tvd(subchannel, dt)
        
        # 5. 应用边界条件
        self.apply_boundary_conditions(subchannel)
        
    def solve_pressure_correction_rhie_chow(self, subchannel):
        """使用Rhie-Chow插值求解压力修正方程"""
        n_axial = len(subchannel.alpha)
        A = np.zeros((n_axial, n_axial))
        b = np.zeros(n_axial)
        
        # 设置边界条件
        A[0, 0] = 1.0
        A[-1, -1] = 1.0
        
        # 构建内部节点方程
        for i in range(1, n_axial-1):
            # 获取物性
            props_g = self.physics.calculate_properties(subchannel.pressure[i], PhaseType.VAPOR)
            props_l = self.physics.calculate_properties(subchannel.pressure[i], PhaseType.LIQUID)
            
            # Rhie-Chow插值系数
            alpha = np.clip(subchannel.alpha[i], 1e-6, 1-1e-6)
            rho_g = max(props_g['density'], 1e-6)
            rho_l = max(props_l['density'], 1e-6)
            
            # 计算动量插值系数
            ap_g = rho_g / self.dt + 2 * props_g['viscosity'] / (self.grid.dz**2)
            ap_l = rho_l / self.dt + 2 * props_l['viscosity'] / (self.grid.dz**2)
            
            # 组装系数矩阵
            coeff = alpha/ap_g + (1-alpha)/ap_l
            A[i, i-1] = -coeff
            A[i, i] = 2 * coeff
            A[i, i+1] = -coeff
            
            # 计算源项（质量不平衡）
            v_g = subchannel.velocity['vapor']['axial'][i]
            v_l = subchannel.velocity['liquid']['axial'][i]
            
            flux_g = alpha * rho_g * v_g
            flux_l = (1-alpha) * rho_l * v_l
            
            b[i] = -(flux_g + flux_l)
        
        # 求解方程组
        try:
            # 添加对角优势
            for i in range(n_axial):
                A[i,i] *= 1.1
            
            # 使用直接求解器
            correction = np.linalg.solve(A, b)
            
            # 限制修正幅度
            max_correction = 0.1 * np.mean(subchannel.pressure)
            correction = np.clip(correction, -max_correction, max_correction)
            
            return correction
            
        except np.linalg.LinAlgError:
            print("警告：压力修正方程求解失败，使用简单梯度估计")
            return np.zeros(n_axial)
            
    def update_void_fraction_tvd(self, subchannel, dt):
        """使用TVD格式更新空泡份额"""
        n = len(subchannel.alpha)
        alpha_new = np.zeros(n)
        
        # 获取物性和漂移速度
        props_g = self.physics.calculate_properties(subchannel.pressure[0], PhaseType.VAPOR)
        props_l = self.physics.calculate_properties(subchannel.pressure[0], PhaseType.LIQUID)
        v_gj, C0 = self.physics.calculate_drift_velocity(props_g['density'], props_l['density'])
        
        # 计算TVD通量限制器
        for i in range(1, n-1):
            # 计算梯度
            r_forward = (subchannel.alpha[i+1] - subchannel.alpha[i]) / \
                       max(subchannel.alpha[i] - subchannel.alpha[i-1], 1e-10)
            
            # 使用van Leer限制器
            phi = (r_forward + abs(r_forward)) / (1 + abs(r_forward))
            
            # 计算有限体积通量
            v_g = subchannel.velocity['vapor']['axial'][i]
            j_drift = subchannel.alpha[i] * (1-subchannel.alpha[i]) * v_gj
            
            # 更新空泡份额
            flux = v_g * subchannel.alpha[i] + j_drift
            alpha_new[i] = subchannel.alpha[i] - dt/self.grid.dz * \
                          (flux * phi - flux * phi)
        
        # 边界条件
        alpha_new[0] = subchannel.alpha[0]
        alpha_new[-1] = alpha_new[-2]
        
        # 限制在物理范围内
        subchannel.alpha = np.clip(alpha_new, 0.0, 1.0)
        
    def apply_explicit_corrections(self, subchannel, dt):
        """应用显式修正项"""
        n = len(subchannel.alpha)
        
        for i in range(1, n-1):
            # 计算扩散项
            for phase in ['vapor', 'liquid']:
                props = self.physics.calculate_properties(
                    subchannel.pressure[i], 
                    PhaseType.VAPOR if phase == 'vapor' else PhaseType.LIQUID
                )
                
                # 计算粘性扩散
                d2v_dz2 = (subchannel.velocity[phase]['axial'][i+1] - 
                          2*subchannel.velocity[phase]['axial'][i] + 
                          subchannel.velocity[phase]['axial'][i-1]) / (self.grid.dz**2)
                
                diffusion = props['viscosity'] / props['density'] * d2v_dz2
                
                # 应用修正
                subchannel.velocity[phase]['axial'][i] += dt * diffusion 